Jednym z najczęstszych równań, od których rozpoczyna się przygodę z metodą elementów skończonych jest równanie transportu ciepła. Naszą niewiadomą funkcją \( u \) jest tutaj funkcja przypisująca wartość temperatury różnym punktom obszaru, w którym rozwiązujemy nasze zadanie. Innymi słowy \( u(x,y) \) reprezentuje wartość temperatury w punkcje \( (x,y) \) obszaru naszej dziedziny \( \Omega \).
Poszukiwana funkcja rozkładu temperatury określona jest w sposób następujący: \( \Omega \ni (x,y) \rightarrow u(x,y)\in R \).

Sformułowanie silne dla problemu transportu ciepła wygląda w sposób następujący
\( -\Delta u=f \)
W równaniu tym \( -\Delta u \) oznacza tak zwany Laplasjan, czyli sumę pochodnych cząstkowych drugiego stopnia, oraz dla uproszczenia zapisu opuściliśmy zmienne \( (x,y) \).
\( \Delta u (x,y) = \frac{\partial^2 u(x,y) } { \partial x^2} + \frac{\partial^2 u(x,y) }{\partial y^2 } \)
natomiast \( f(x,y) \) oznacza funkcję reprezentującą źródła ciepła. W tych fragmentach dziedziny, w których funkcja \( f(x,y) \) jest dodatnia, znajduje się "ogrzewanie", natomiast w tych fragmentach dziedziny, w których funkcja \( f(x,y) \) jest ujemna, znajduje się "chłodzenie".

Sformułowanie słabe uzyskuje się w sposób następujący. Całkujemy i przemnażamy nasze równanie przez wybrane funkcje \( v(x,y) \) zwane funkcjami testującymi, których używać będziemy do uśredniania naszego równania w obszarze, na którym funkcje te są określone
\( -\int_{\Omega} \Delta u (x,y) v(x,y) dxdy = \int_{\Omega} f(x,y) v(x,y) dxdy \)

W podobny sposób jak miało to miejsce w przypadku problemu projekcji bitmapy, każdy wybór funkcji testującej \( v(x,y) \) służącej do uśredniania naszego problemu w obszarze, w którym funkcja testująca jest określona, daje nam jedno równanie. Różne wybory funkcji testującej \( v(x,y) \) dadzą nam więc różne równania .

Jeśli rozpiszemy nasz wzór na pochodne cząstkowe, dostajemy dwie osobne całki
\( \int_{\Omega} \Delta u(x,y)v(x,y) dxdy= \int_{\Omega} \left( \frac{\partial^2 u(x,y) }{\partial x^2 } + \frac{\partial^2 u(x,y) }{\partial y^2 } \right) v(x,y) dxdy \)

Nasz problem, który chcemy rozwiązać wygląda więc tak: Szukamy rozkładu temperatury, funkcji \( \Omega \ni (x,y) \leftarrow u(x,y) \in R \) takiej że
\( \int_{\Omega} \frac{\partial^2 u(x,y)}{\partial x^2} v(x,y) dxdy + \int_{\Omega} \frac{\partial^2 u(x,y)}{\partial y^2} v(x,y) dxdy = \int_{\Omega} f(x,y) v(x,y) dxdy \)
dla dowolnych funkcji testujących \( \Omega \ni (x,y) \rightarrow v(x,y) \in R \) (oczywiście nasz rozkład temperatury oraz funkcja testująca muszą być takie, żeby dało się policzyć całki, jak wiemy, nie ze wszystkich funkcji da się policzyć całki, niektóre funkcje dają całki, które mogą być nieskończone),
\( a(u,v)=l(v) \\ a(u,v)=\int_{\Omega } \frac{\partial^2 u(x,y) }{\partial x^2 } v(x,y) dxdy + \int_{\Omega } \frac{\partial^2 u(x,y) }{\partial y^2 } v(x,y) dxdy \\ l(v)=\int_{\Omega} f(x,y) v(x,y) dxdy \)

Zauważmy teraz, że gdybyśmy chcieli zastosować teraz nasze funkcje B-spline wprowadzone w poprzednim rozdziale do przybliżenia pola temperatury,
\( u(x,y) = \sum_{i=1,j=1}^{N_x,N_y} u_{i,j} B^x_{i}(x) B^y_{j}(y) \)
wówczas musielibyśmy być w stanie policzyć pochodne drugiego stopnia z funkcji B-spline.

Żeby pochodne drugiego stopnia były możliwe do policzenia, funkcje bazowe użyte do przybliżenia pola temperatury muszą być smukłe. B-spline'y świetnie nadają się do liczenia z nich pochodnych. Jednakże, żeby pochodne te nie wyszły równe zero (żeby nasz układ równań, który buduje się podobnie jak w przypadku przybliżania bitmapy nie wyszedł równy zero), B-spline'y stosowane tutaj muszą być funkcjami wyższego rzędu (co najmniej klasy \( C^2 \), czyli muszą być co najmniej B-spline'ami trzeciego stopnia).
Jeśli lubimy wielomiany wyższego rzędu, możemy w tym momencie wstawić całki naszego równania do macierzy problemu, czyli podmienić macierz problemu projekcji bitmapy naszą macierzą.

Teraz, w miejsce funkcji \( u \) wstawiamy przybliżenie za pomocą funkcji B-spline \( u = \sum_{i=1,j=1}^{N_x,N_y} u_{i,j} B^x_{i} B^y_{j} \), oraz w miejsce funkcji testujących również wstawiamy funkcje B-spline \( B^x_{k} B^y_{l}, k=1,...,N_x; j=1,...,N_y \). Analogicznie jak w przypadku projekcji bitmapy dostajemy układ równań
\( \begin{bmatrix} a({B^x_{1,p}B^y_{1,p}},B^x_{1,p}B^y_{1,p}) & a({B^x_{1,p}B^y_{1,p}},{B^x_{2,p}B^y_{1,p}}) & \cdots & a({B^x_{1,p}B^y_{1,p}},{B^x_{N_x,p}B^y_{N_y,p}}) \\ a({B^x_{2,p}B^y_{1,p}},{B^x_{1,p}B^y_{1,p}}) & a({B^x_{2,p}B^y_{1,p}}{B^x_{2,p}B^y_{1,p}}) & \cdots & a({B^x_{2,p}B^y_{1,p}},{B^x_{N_x,p}B^y_{N_y,p}}) \\ \vdots & \vdots & \vdots & \vdots \\ a({B^x_{N_x,p}B^y_{N_y,p}},{B^x_{1,p}B^y_{1,p}}) & a({B^x_{N_x,p}B^y_{N_y,p}},{B^x_{2,p}B^y_{1,p}}) & \cdots & a({B^x_{N_x,p}B^y_{N_y,p}},{B^x_{N_x,p}B^y_{N_y,p}}) \\ \end{bmatrix} \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ \vdots \\ u_{N_x,N_y}\\ \end{bmatrix} \\ = \begin{bmatrix} l( B^x_1(x)*B^y_1(y)) \\ l(B^x_1(x)*B^y_2(y)) \\ \vdots \\ l(B^x_{N_x}(x)*B^y_{N_y}(y)) \\ \end{bmatrix} \)
gdzie po uwzględnieniu wzorów na \( a \) i \( l \) mamy
\( \begin{bmatrix} \int{\Delta( B^x_{1,p}B^y_{1,p }) }{B^x_{1,p}B^y_{1,p} } & \int{\Delta( B^x_{1,p}B^y_{1,p} ) }{B^x_{2,p}B^y_{1,p} } & \cdots & \int{\Delta( B^x_{1,p}B^y_{1,p}) }{B^x_{N_x,p}B^y_{N_y,p} } \\ \int{ \Delta( B^x_{2,p}B^y_{1,p}) }{ B^x_{1,p}B^y_{1,p} } & \int{ \Delta(B^x_{2,p}B^y_{1,p}) }{ B^x_{2,p}B^y_{1,p} } & \cdots & \int{ \Delta( B^x_{2,p}B^y_{1,p} ) }{ B^x_{N_x,p}B^y_{N_y,p} } \\ \vdots & \vdots & \vdots & \vdots \\ \int{ \Delta( B^x_{N_x,p}B^y_{N_y,p} ) }{B^x_{1,p}B^y_{1,p} } & \int{ \Delta( B^x_{N_x,p }B^y_{N_y,p}) }{B^x_{2,p}B^y_{1,p} } & \cdots & \int{ \Delta( B^x_{N_x,p }B^y_{N_y,p} ) }{B^x_{N_x,p}B^y_{N_y,p} } \\ \end{bmatrix} \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ \vdots \\ u_{N_x,N_y}\\\end{bmatrix} \\ = \begin{bmatrix} \int f(x,y) B^x_1(x)*B^y_1(y) dx \\ \int f(x,y) B^x_1(x)*B^y_2(y) dx \\ \vdots \\ \int f(x,y) B^x_{N_x}(x)*B^y_{N_y}(y) dx \\ \end{bmatrix} \)

Oczywiście musimy pamiętać, że
\( \int_{\Omega} \Delta u(x,y)v(x,y) dxdy= \int_{\Omega} \left( \frac{\partial^2 u(x,y) }{ \partial x^2 } + \frac{\partial^2 u(x,y) }{ \partial y^2 } \right) v(x,y) dxdy \)

Teraz, gdybyśmy wygenerowali stosowne całki i uruchomili jeden z algorytmów faktoryzacji macierzy, opisany w rozdziale piątym, okazałoby się, że macierz jest osobliwa i układu równań nie da się rozwiązać! Dlaczego tak by się stało? Ponieważ kiedy rozwiązujemy problem transportu ciepła, i określamy równanie spełnione na obszarze \( \Omega \) musimy sprecyzować tak zwane warunki brzegowe, czyli opisać jak temperatura zachowuje się na brzegu obszaru.

Zawsze kiedy chcemy sformułować nasz problem inżynieryjny za pomocą równań różniczkowych cząstkowych, na przykład problem transportu ciepła, na jakimś obszarze \( \Omega \), musimy po pierwsze podać równanie różniczkowe cząstkowe, które opisuje modelowane zjawisko wewnątrz obszaru \( \Omega \), oraz dodatkowo musimy opisać tak zwane warunki brzegowe, czyli podać pewne dodatkowe równania, które opisują zachowanie naszego zjawiska na brzegu obszaru.

Wyobraźmy sobie, że chcemy obliczyć temperaturę we wnętrzu pokoju. Musimy oczywiście sprecyzować jakie równanie różniczkowe cząstkowe opisuje zjawisko fizyczne przepływu ciepła wewnątrz pokoju. Jest to oczywiście znane nam już równanie (sformułowanie silne) z Laplasjanem \( \Delta u(x,y)=f(x,y) \). Jeśli podamy samo równanie spełnione w środku pokoju, ale nie podamy tego co dzieje się na ścianach pokoju, mamy potencjalnie nieskończenie wiele możliwych rozwiązań. Naszego układ równań nie da się więc rozwiązać bez sprecyzowania tego co dzieje się na brzegu. Poprawne określenie warunków brzegowych daje nam zazwyczaj jedno możliwe rozwiązanie.

Nasz obszar \( \Omega \) może więc reprezentować wnętrze pokoju. Dla uproszczenia przyjmijmy że pracujemy w dwóch wymiarach (że nasz pokój jest płaski dwuwymiarowy). Załóżmy, że pokój ten ma kształt litery L, oraz jego wewnętrzne ściany przylegają do chłodnicy, na brzegu której temperatura wynosi zero, a pozostałe ściany mają okna, przez które dociera do nas ciepło zewnętrzne. Wówczas temperatura na wewnętrznych ścianach pokoju wyniesie zero. Na części brzegu obszaru (na ścianach wewnętrznych, na których zmierzyliśmy temperaturę) wprowadzamy tak zwany warunek brzegowy Dirichleta. Sciany te (część brzegu \( \partial \Omega \)) oznaczamy symbole \( \Gamma_D \)i podajemy wartości temperatury na tym fragmencie brzegu (na ścianach obszaru)
\( u(x,y) = T(x,y) \textrm{ dla }(x,y)\in \Gamma_D \)
gdzie \( T(x,y) \) to jest temperatura zmierzona termometrem lub kamerą termowizyjną. Sytuacja ta jest zilustrowana na rysunku Rys. 1.

Na pewnym fragmencie brzegu pokoju mogą znajdować się również okna, a za oknami może być gorąco (jeśli jest lato) lub zimno (jeśli jest zima) i temperatura wokół szyby może cały czas rosnąć (w lecie) lub maleć (w zimie). Wówczas za pomocą termometru i zegarka (lub magicznej kamery termowizyjnej) możemy zmierzyć prędkość z jaką zmienia się temperatura obok szyby.
W takim przypadku na części brzegu obszaru (na oknach, na których zmierzyliśmy prędkość zmian temperatury) wprowadzamy tak zwany warunek brzegowy Neumanna. Ten fragment brzegu (czyli okna) oznaczamy symbolem \( \Gamma_N \) i określamy prędkość zmian temperatury
\( \frac{\partial u(x,y)}{\partial n} = g(x,y) \textrm{ dla }(x,y)\in \Gamma_N \)
gdzie \( g(x,y) \) to strumień ciepła na brzegu obszaru.
Co oznacza tutaj \( \frac{\partial u }{\partial n } \)? Oznacza prędkość zmian temperatury w kierunku brzegu. Z matematycznego punktu widzenia \( \frac{\partial u }{\partial n }=\frac{\partial u }{\partial x }n_1+\frac{\partial u }{\partial y }n_2 \) gdzie \( n=(n_1,n_2) \) to wektor prostopadły do brzegu obszaru \( \Omega \), skierowany na zewnątrz, czyli np. jeśli brzeg to ściana prostopadła do osi \( x \) to \( n=(1,0) \), jeśli brzeg to dolna ściana pokoju, to \( n=(0,-1) \) itd. Przykładowe wektory \( n=(n_1,n_2) \) przedstawiono na rysunku Rys. 1.

W jaki sposób zmodyfikować nasz uklad równań tak żeby uwzględnić w nim zerową temperaturę na części brzegu \( \Gamma_D \) oraz prędkość przepływu ciepła na części brzegu \( \Gamma_N \).
W tym celu musimy sprecyzować jak wyglądać będzie nasza grupa elementów oraz nasze funkcje bazowe B-spline rozpięte na elementach, podając wektor węzłów wzdłuż osi \( x \) oraz wektor węzłów wzdłuż osi \( y \). Najprostszym rozwiązaniem będzie podzielenie naszego obszaru na cztery podobszary kwadratowe (cztery grupy elementów) i sklejenie ich ze soba (powtórzenie węzłów na interfejsie sklejenia tak żeby funkcje bazowe w tym miejscu miały ciągłość C0). Odsyłamy tutaj do rozdziału trzeciego, w którym opisany jest sposób definiowania funkcji bazowych za pomocą wektorów węzłów i wielomianów B-spline.
Wzdłuż osi \( x \) wprowadzamy wektor węzłów [0 0 0 1 2 2 3 4 4 4], podobnie wzdłuż osi y wprowadzamy wektor węzłów [0 0 0 1 2 2 3 4 4 4].

Powtórzyliśmy środkowy węzeł tak żeby zredukować ciągłość funkcji B-spline na środku obszaru. Pozwoli nam to wymusić poprawnie warunek brzegowy Dirichleta na jednej ćwiartce obszaru. Zauważmy że funkcje B-spline są smukłe i rozpinają się na wiele elementów, nie da się więc bez redukcji ciągłości (bez powtarzania węzłów w wektorach węzłów) wymuszać wartości funkcji B-spline w jakimś punkcie.
Wynikające z nich funkcje bazowe przedstawione są na rysunku Rys. 2. Uzyskaliśmy dwie bazy jednowymiarowych funkcji B-spline
\( \{ B_{i,2}(x) \}_{i=1,...,5 } \) oraz \( \{B_{j,2}(y)\}_{j=1,...,5 } \). Następnie utworzymy z nich iloczyny tensorowe \( B_{i,j;2}(x,y)=B_{i,2}(x)B_{j,2}(y),i,j=1,...,5 \).

Zauważmy, że nasz obszar \( \Omega \) rozpięty jest na kwadracie \( [-1,1]^2 \). Natomiast nasze 5*5=25 funkcji bazowych rozpiętych jest na kwadracie \( [0,2]^2 \) wynikającym z definicji wektora węzłów [0 0 1 2 2 3 4 4].

Musimy więc "przemapować" nasz obszar fizyczny \( [-1,1]^2 \) na grupę elementów, na której zdefiniowane są nasze funkcje B-spline
\( [-1,1]^2 \ni (x,y) \rightarrow \psi(x,y)= \left(2(x+1),2(y+1)\right) \)
Rozszerzamy nasze funkcje na obszar \( [-1,1]^2 \) zakładając, że poza naszym obszarem są one równe zero, oraz przeprowadzamy zmianę zmiennych w całce z bloku \( [-1,1]^2 \) na blok \( [0,4]^2 \)

Wówczas wszystkie całki przemapowujemy na nasz obszar, na którym rozpinamy całki
\( a(u,v)=\int_{\Omega} \frac{\partial^2 u(x,y) }{ \partial x^2 } v(x,y) dxdy + \int_{\Omega} \frac{\partial^2 u(x,y) }{ \partial y^2 } v(x,y) dxdy= \int_{[-1,1]^2} \frac{\partial^2 u(x,y) }{ \partial x^2 } v(x,y) dxdy + \int_{[-1,1]^2 } \frac{\partial^2 u(x,y) }{ \partial y^2 } v(x,y) dxdy= \\ \int_{[0,4]^2} \frac{\partial^2 u(x,y) }{\partial x^2 } v(x,y) |Jac(\psi^{-1})| dxdy + \int_{[0,4]^2} \frac{\partial^2 u(x,y)}{ \partial y^2 } v(x,y) |Jac(\psi^{-1})| dxdy \\ l(v)=\int_{\Omega} f(x,y) v(x,y) dxdy = \\= \int_{[-1,1]^2} f(x,y) v(x,y) dxdy = \int_{[0,4]^2} f(x,y) v(x,y) |Jac(\psi^{-1})| dxdy \)

W powyższym wzorach na całki umieściliśmy wartość bezwzględną z Jakobianu transformacji zmiany zmiennych z naszego obszaru \( \Omega=[-1,1]^2 \) na obszar \( [0,4]^2 \), na którym rozpięliśmy wektory węzłów. Jest ona obliczana tak. Znamy odwzorowanie \( \psi:\Omega \rightarrow [0,4]^2 \) dane \( \psi(x,y)=\left(2(x+1),2(y+1)\right) \). Oznacza to, że nowe zmienne w obszarze \( [0,4]^2 \) (które oznaczmy na chwile \( (\xi,\eta) \) dane są wzorami \( (\xi,\eta)= \left(2(x+1),2(y+1)\right) \). Jeśli chciałbym obliczyć zmienne \( (x,y) \) w obszarze \( [-1,1]^2 \) to muszę po prostu przeliczyć
\( 2(x+1)=\xi \) czyli \( x=\xi/2-1 \), oraz \( 2(y+1)=\eta \), czyli \( y=\eta/2-1 \).
Tak więc \( \psi^{-1}(\xi,\eta)=(\xi/2-1,\eta/2-1) \). Dla wygody zmieńmy zmienne z powrotem na symbole \( (x,y) \) czyli \( \psi^{-1}(\xi,\eta)=(x/2-1,y/2-1) \). Teraz, musimy policzyć pochodne cząstkowe \( \frac{\partial(x/2-1)}{\partial x }=1/2 \), \( \frac{\partial(x/2-1) }{\partial y }=0 \),
\( \frac{\partial(y/2-1) }{\partial x }=0 \),
oraz \( \frac{\partial(y/2-1) }{\partial y }=1/2 \). Nasz Jakobian w całce to wyznacznik z macierzy pochodnych cząstkowych
\( jac(\psi^{-1})=|\begin{bmatrix} 1/2 & 0 \\ 0 & 1/2 \end{bmatrix}| = |1/2*1/2|=|1/4|=1/4 \)
czyli nasz układ równań
\( a(u,v)=\int_{[0,4]^2} \frac{\partial^2 u(x,y) }{ \partial x^2 } v(x,y) 1/4 dxdy + \int_{ [0,4]^2 } \frac{ \partial^2 u(x,y) }{\partial y^2 } v(x,y) 1/4 dxdy \\ l(v)=\int_{[0,4]^2 } f(x,y) v(x,y) 1/4 dxdy \)

Teraz, dla funkcji B-spline rozpiętych na brzegu Dirichleta obszaru możemy wymusić warunek Dirichleta w ten sposób, że zerujemy stosowne wiersze i prawe strony układu równań, oraz na przekątnej wpisujemy jedynki.
Innymi słowy, dla stosownych wierszy odpowiadających B-spline'om \( i,j \) położonym na brzegu Dirichleta
\( \int{\Delta(B^x_{i,p }B^y_{j,p } ) }{B^x_{1,p } B^y_{1,p } }u_{1,1 } + \cdots + \int{ \Delta(B^x_{i,p }B^y_{j,p}) }{ B^x_{2,p }B^y_{1,p} }u_{i,j }+ \cdots +\int{ \Delta( B^x_{i,p }B^y_{j,p }) }{B^x_{N_x,p }B^y_{N_y,p } }u_{N_x,N_y } = \\= \int f(x,y)B^x_{ i,p }B^y_{j,p } dxdy \\ \)
zamieniamy je na
\( 1*u_{i,j }=0 \)

Jak natomiast wyrazić warunek Neumana? Nie jest to niestety łatwo wykonalne, bez całkowania przez części.
Musimy więc zastosować wzór na całkowanie przez części. Ponieważ nasz problem jest dwuwymiarowy, wzór ten jest stosunkowo skomplikowany. Nie będziemy go tutaj wyprowadzać, ponieważ wymagało by to przeprowadzenia kursu z analizy matematycznej.
\( \int_{\Omega} \Delta u (x,y) v(x,y) dxdy =- \int_{\Omega} \nabla u (x,y) \cdot \nabla v(x,y) dxdy + \int_{\partial \Omega} \frac{\partial u }{\partial n } v dS \)
gdzie \( \nabla u(x,y) = (\frac{\partial u}{\partial x }, \frac{\partial u }{\partial y }) \) oznacza gradient \( u \) czyli wektor pochodnych cząstkowych funkcji \( u \), podobnie \( \nabla v(x,y) = (\frac{\partial v}{\partial x }, \frac{\partial v }{\partial y }) \)
oznacza gradient \( v \) czyli wektor pochodnych cząstkowych funkcji \( v \), natomiast \( \cdot \) oznacza iloczyn skalarny ich wektorów, czyli
\( (\frac{\partial u }{\partial x }, \frac{\partial u }{\partial y }) \cdot (\frac{\partial v }{\partial x }, \frac{\partial v }{\partial y }) = \frac{\partial u }{\partial x }\frac{\partial v }{\partial x } + \frac{\partial u }{\partial y } \frac{\partial v }{\partial y } \).

Tak więc całe nasze równanie
\( -\int_{\Omega} \frac{\partial^2 u(x,y)}{\partial x^2 } v(x,y) dxdy - \int_{\Omega} \frac{\partial^2 u(x,y)}{\partial y^2} v(x,y) dxdy = \int_{\Omega} f(x,y) v(x,y) dxdy \)
po całkowaniu przez części zmienia się w
\( \int_{\Omega} \left(\frac{\partial u }{\partial x }\frac{\partial v }{\partial x } + \frac{\partial u }{\partial y } \frac{\partial v }{\partial y } \right) dxdy= \int_{\partial \Omega } \frac{\partial u }{\partial n } v dS +\int_{\Omega } f(x,y) v(x,y) dxdy \)
Co oznacza tutaj \( \frac{\partial u }{\partial n } \)? Oznacza prędkość zmian temperatury w kierunku brzegu. Z matematycznego punktu widzenia \( \frac{\partial u }{\partial n}=\frac{\partial u }{\partial x }n_1+\frac{\partial u }{\partial y }n_2 \) gdzie \( n=(n_1,n_2) \) to wektor prostopadły do brzegu obszaru \( \Omega \), czyli np. jeśli brzeg to pionowa ściana prostopadła do osi \( x \) to \( n=(1,0) \), jeśli brzeg to podłoga pokoju, to \( n=(0,-1) \) itd.
Dziedzina, na której rozwiązujemy nasz problem \( \Omega = [-1,1] \times [-1,1] \ [-1,0] \times [-1,0] \) ma kształt odwróconej litery L.

Na wewnętrznej części brzegu obszaru, którą oznaczamy \( \Gamma_D=\{ (x,y): x=0, y\in[-1,0] \textrm{ lub } x \in [-1,0], y=0 \} \) przyjmujemy zerową temperaturę wyrażoną warunkiem brzegowym Dirichleta
\( u(x,y) = 0 \textrm{ dla }(x,y)\in \Gamma_D \)
Na pozostałej części brzegu obszaru, którą oznaczamy \( \Gamma_N=\Gamma \ \Gamma_D \) wprowadzamy tak zwany warunek brzegowy Neumanna.
Ten fragment brzegu (czyli okna) oznaczamy symbolem \( \Gamma_N \) i określamy prękość zmian temperatury
\( \frac{\partial u(x,y)}{\partial n} = g(x,y) \textrm{ dla }(x,y)\in \Gamma_N \)
gdzie \( g(x,y) \) to zadana prędkość zmian temperatury na brzegu.

Podsumowując, brzeg naszego obszaru \( \partial \Omega \) podzieliliśmy na dwa fragmenty \( \partial \Omega = \Omega_D \cup \Omega_N \) i na pierwszym zmierzyliśmy temperaturę, na drugim zmierzyliśmy prędkość zmian temperatury.
Możemy więc fragment naszego sformułowania słabego z całką po brzegu rozbić na te trzy fragmenty brzegu
\( \int_{\partial \Omega } \frac{\partial u(x,y) }{\partial n } v dS = \int_{\partial \Omega_D } \frac{ \partial u(x,y) }{\partial n } v(x,y) dS + \\ \int_{\partial \Omega_N } \frac{\partial u(x,y) }{\partial n } v(x,y) dS \)

Co zrobić dalej z tymi trzema całkami rozpiętymi na trzech fragmentach brzegu?
Całka z warunkiem brzegowym Neumanna jest prosta, ponieważ z naszego warunku brzegowego Neumanna wynika, że
\( \frac{\partial u(x,y) }{\partial n }=g(x,y) \), gdzie \( g(x,y) \) to dana prędkość zmian temperatury, czyli \( \int_{\partial \Omega_N } \frac{ \partial u(x,y) }{\partial n } v(x,y) dS=\int_{\partial \Omega_N } g(x,y)v(x,y) dS \). Dlaczego musimy zamienić w tej całce pochodną w kierunku normalnym na zadaną funkcję \( g \)? Ponieważ chcemy zadać ten warunek brzegowy, chcemy zamodelować w naszej symulacji sytuację kiedy ciepło naprawdę w taki sposób oddziaływuje przez okno na wnętrze naszego pokoju.

Jeśli chodzi o warunek brzegowy Dirichleta, to wymaga on pewnych dodatkowych zabiegów.
Jeśli \( T(x,y)=0 \) w warunku Dirichleta, czyli \( u(x,y) = 0 \textrm{ dla }(x,y)\in \Gamma_D \) wówczas sprawa jest prosta. Wiemy, że na tym brzegu temperatura wynosi zero, nie musimy więc obliczać średnich wartości w tym miejscu brzegu, nie musimy wybierać takich funkcji testujących \( v(x,y) \), które nie są równe zero na brzegu Dirichleta. Innymi słowy mogę przyjąć, że funkcje testujące \( v(x,y)=0 \) dla \( (x,y) \in \Gamma_D \). Innymi słowy wówczas całka
\( \int_{\partial \Omega_D } \frac{\partial u(x,y) }{\partial n } v(x,y) dS=0 \).
Nasza prawa strona
\( l(v) = \int_{\partial \Omega_D } \frac{\partial u(x,y) }{\partial n } v(x,y) dS + \int_{\partial \Omega_N } \frac{\partial u(x,y) }{\partial n } v(x,y) dS \\ +\int_{\Omega } f(x,y) v(x,y) dxdy = \int_{\partial \Omega_N } g(x,y) v(x,y) dS +\int_{\Omega } f(x,y) v(x,y) dxdy \)
Nasz problem wygląda więc tak
\( a(u,v)=l(v) \\ a(u,v) = \int_{\Omega} \left(\frac{\partial u }{\partial x }\frac{\partial v }{\partial x } + \frac{\partial u }{\partial y } \frac{\partial v }{\partial y } \right)v(x,y) dxdy \\ l(v) = \int_{\partial \Omega_N } g(x,y) v(x,y) dS \nonumber +\int_{\Omega } f(x,y) v(x,y) dxdy \)

Zauważmy, że skoro całkowaliśmy przez części, wówczas możemy pozwolić sobie na redukcję stopnia funkcji B-spline (teraz, nasz wzór zawiera jedynie pochodne pierwszego rzędu, i wystarczą nam funkcje klasy \( C^1 \)).
Ponadto, im wyższy stopień B-spline'ów, tym bardziej kosztowne jest jego przetwarzanie. Samo wygenerowanie wzoru na B-spline wyższego rzędu wymaga więcej operacji, podobnie obliczanie wartości takiego B-spline'a. Ponadto układ równań wygenerowany dla B-spline'ów wyższego rzędu będzie gęstszy, będzie miał mniej zer, i przez to koszt faktoryzacji (koszt rozwiązania takiego układu równań) będzie większy.